*************************************************************;
***figure_7.sas                                           ***;
*************************************************************;
***This program creates data for Figure 7 for the Davis,  ***;
***Grim, Haltiwanger, and Streitwieser 2013 RESTAT paper. ***;
*************************************************************;
***Output:                                                ***;
*** ~/csv/figure7_1a.csv                                  ***;
*** ~/csv/figure7_2.csv                                   ***;
***Dataset for later use:                                 ***;
*** elec.styr_elas                                        ***;
*************************************************************;

options ls=80 obs=max;

****MACROS****;

**Macros for top panels of Figure 7;

*Create weighted mean elasticities;
	
%macro d1(ds, wt, var);

proc sort data=&ds;
	by fipsst year;

proc summary data=&ds noprint;
	by fipsst year;
	var &var;
	weight wt&wt;
	output out=mean_&var._&wt
	mean(&var) = mean_&var._&wt;

%mend d1;

%macro w(year);

data w&year (keep = fipsst db_adj_&year ddb_adj_&year efr_adj_&year
			fac_i_adj_&year fr_adj_&year inter_adj_&year
			rank_s_adj_&year rank_s_short_adj_&year seas_adj_&year);
	set naruc;
	if year = &year;
	db_adj_&year = db_adj;
	ddb_adj_&year = ddb_adj;
	efr_adj_&year = efr_adj;
	fac_i_adj_&year = fac_i_adj;
	fr_adj_&year = fr_adj;
	inter_adj_&year = inter_adj;
	rank_s_adj_&year = rank_s_adj;
	rank_s_short_adj_&year = rank_s_short_adj;
	seas_adj_&year = seas_adj;
	
proc sort data=w&year;
	by fipsst;

%mend w;

%macro d(year);

data e&year (keep = fipsst postalst elas_: sum_pe_&year wtpe100_&year wttvs100_&year);
	set e;
	if year = &year;
	elas_nf_pe_&year = mean_elas_nf_pe;
	elas_pe_&year = mean_elas_pe;
	elas_nf_tvs_&year = mean_elas_nf_tvs;
	elas_tvs_&year = mean_elas_tvs;
	sum_pe_&year = sum_pe;
	wtpe100_&year = wtpe100;
	wttvs100_&year = wttvs100;

proc sort data=e&year;
	by fipsst;

%mend d;

*Macro to grab US total PE for a particular year;

%macro g(year);

data us&year (keep = us&year);
	set ustot;
	if year = &year;
	us&year = us_pe;

%mend g;

%macro r(var);

data eboth_&var;
	set eboth;

proc sort data=eboth_&var;
	by descending d&var;

data eboth_&var;
	retain pe_tot;
	set eboth_&var;

	rank_&var = _N_;
	
	if _N_ = 1 then pe_tot = sum_pe_1967;
	else pe_tot = pe_tot + sum_pe_1967;

	pct_tot_top = pe_tot / us1967 * 100;	

	if pct_tot_top > 25 then flag_top1 = 1;
	else flag_top1 = 0;

data topn;
	set eboth_&var;
	if flag_top1 = 1;

proc sort data=topn;
	by rank_&var;

data topn (keep = topn_cutoff top_cut_percent);
	set topn;
	if _N_ = 1;
	topn_cutoff = rank_&var;
	top_cut_percent = pct_tot_top;

data eboth_&var;
	if _N_ = 1 then set topn;
	set eboth_&var;

proc print data=eboth_&var;
	var postalst rank_&var topn_cutoff sum_pe_1967 us1967 pct_tot_top flag_top1 d&var &var._1967 &var._1982 ;
	title1 "Check eboth:  &var";

proc sort data=eboth_&var;
	by descending rank_&var;

data eboth_&var;
	retain pe_tot2;
	set eboth_&var (drop = pe_tot flag_top1 pct_tot_top);

	if _N_ = 1 then pe_tot2 = sum_pe_1967;
	else pe_tot2 = pe_tot2 + sum_pe_1967;

	pct_tot_bottom = pe_tot2 / us1967 * 100;	

	if pct_tot_bottom > 25 then flag_bottom1 = 1;
	else flag_bottom1 = 0;

data bottomn;
	set eboth_&var;
	if flag_bottom1 = 1;

proc sort data=bottomn;
	by descending rank_&var;

data bottomn (keep = bottomn_cutoff bottom_cut_percent);
	set bottomn;
	if _N_ = 1;
	bottomn_cutoff = rank_&var;
	bottom_cut_percent = pct_tot_bottom;

data eboth_&var;
	if _N_ = 1 then set bottomn;
	set eboth_&var;

proc print data=eboth_&var;
	var postalst rank_&var bottomn_cutoff sum_pe_1967 us1967 pct_tot_bottom flag_bottom1 d&var &var._1967 
&var._1982 ;
	title1 "Check eboth:  &var";

data eboth_&var;
	set eboth_&var (drop = pct_tot_bottom flag_bottom1 pe_tot2);
	if rank_&var le topn_cutoff then top_&var = 1;
	else top_&var = 0;
	rename topn_cutoff = top_cut_&var;
	rename top_cut_percent = top_cutp_&var;

	if rank_&var ge bottomn_cutoff then bottom_&var = 1;
	else bottom_&var = 0;
	rename bottomn_cutoff = bottom_cut_&var;
	rename bottom_cut_percent = bottom_cutp_&var;

proc sort data=eboth_&var;
	by postalst;

%mend r;


*Create weighted mean elasticities;
	
%macro rd1(ds, wt, var);

%*Top;

data top&ds;
	set &ds;
	if top_&var._&wt = 1;

proc sort data=top&ds;
	by year;

proc summary data=top&ds noprint;
	by year;
	var &var;
	weight wt&wt;
	output out=top_mean_&var._&wt
	mean(&var) = top_mean_&var._&wt;

%*Bottom;

data bottom&ds;
	set &ds;
	if bottom_&var._&wt = 1;

proc sort data=bottom&ds;
	by year;

proc summary data=bottom&ds noprint;
	by year;
	var &var;
	weight wt&wt;
	output out=bottom_mean_&var._&wt
	mean(&var) = bottom_mean_&var._&wt;

data f&var._&wt;
	merge top_mean_&var._&wt bottom_mean_&var._&wt;
	by year;

proc datasets lib=work;
	delete top&ds bottom&ds top_mean_&var._&wt bottom_mean_&var._&wt;
run;
quit;

%mend rd1;

**Macros - bottom panels of Figure 7;

%macro s(num);

proc summary data=all;
	output out=q&num
	p25(index&num) = p25_&num
	p50(index&num) = p50_&num
	p75(index&num) = p75_&num
	;

data all;
	if _N_ = 1 then set q&num;
	set all;

	if index&num < p25_&num then q&num = 1;
	else if index&num ge p25_&num and index&num < p50_&num then q&num = 2;
	else if index&num ge p50_&num and index&num < p75_&num then q&num = 3;
	else if index&num ge p75_&num then q&num = 4;

	if q&num = 4 then top&num = 1;
	else top&num = 0;

	if q&num = 1 then bottom&num = 1;
	else bottom&num = 0;

%*Get data for csv of top and bottom states by each measure;

data top&num (keep = postalst what);
	length what $20.;
	set all;
	if top&num = 1;

	what = "top&num";

data bottom&num (keep = postalst what);
	length what $20.;
	set all;
	if bottom&num = 1;

	what = "bottom&num";

data both&num;
	set top&num bottom&num;

%mend s;


*Create weighted mean elasticities;
	
%macro d1_2(ds, wt, var);

%*Top;

data top&ds;
	set &ds;
	if top1 = 1;

proc sort data=top&ds;
	by year;

proc summary data=top&ds noprint;
	by year;
	var &var;
	weight wt&wt;
	output out=top_mean_&var._&wt.1
	mean(&var) = top_mean_&var._&wt.1;

%*Bottom;

data bottom&ds;
	set &ds;
	if bottom1 = 1;

proc sort data=bottom&ds;
	by year;

proc summary data=bottom&ds noprint;
	by year;
	var &var;
	weight wt&wt;
	output out=bottom_mean_&var._&wt.1
	mean(&var) = bottom_mean_&var._&wt.1;

data f&var._&wt.1;
	merge top_mean_&var._&wt.1 bottom_mean_&var._&wt.1;
	by year;

proc datasets lib=work;
	delete top&ds bottom&ds top_mean_&var._&wt.1 bottom_mean_&var._&wt.1;
run;
quit;

%mend d1_2;

****START PROGRAM****;

data outurall;
	set elec.outurall;
	pe = exp(logpe);

data outurall_nf;
	set elec.outurall_nf;

%d1(outurall, tvs, elas)
%d1(outurall, pe, elas)
%d1(outurall_nf, tvs, elas_nf)
%d1(outurall_nf, pe, elas_nf)

*Calculate state-year sum tvs and pe (sample-weighted);

proc summary data=outurall;
	by fipsst year;
	weight wt;
	output out=sum
	sum(pe) = sum_pe
	sum(tvs) = sum_tvs
	;

*Create state purchase-weighted mean price (unit value);

proc summary data=outurall vardef=wdf;
        by fipsst year;
        weight wtpe;
        output out=meanp
        mean(lpe_gdp) = mean_lpe_gdp_st;

data styr_elas1;
	merge mean_elas_tvs mean_elas_pe mean_elas_nf_tvs mean_elas_nf_pe sum meanp;
	by fipsst year;

proc print data=styr_elas1;
	title1 "Elasticity Dataset";

*State-year mean elasticities;

data pqem;
	set styr_elas1 (keep = fipsst year mean_elas_nf_tvs mean_elas_tvs
		mean_elas_nf_pe mean_elas_pe sum_pe sum_tvs);

proc sort data=pqem;
	by fipsst year;

*Fuel shares and fossil fuel cost index (simple impute - impute 2);

data seds;
	set eia.seds2008select (keep = fipsst year share_gen: pf1_2 pf2_2 pf3_2 censt postalst state division region); *Select variables from 2008 EIA State Energy Data System data 
available at www.eia.doe.gov with simple imputes for missing values; 
	rename pf1_2 = PF1;
	rename pf2_2 = PF2;
	rename pf3_2 = PF3;

data seds (drop = prev_:);
	set seds;

	*Calculate 1 yr log change PF variables;
	if year ne 1960 then do;
		prev_pf1 = lag(pf1);
		prev_pf2 = lag(pf2);
		prev_pf3 = lag(pf3);

		dpf1 = (log(pf1) - log(prev_pf1));
		dpf2 = (log(pf2) - log(prev_pf2));
		dpf3 = (log(pf3) - log(prev_pf3));
	end;

proc sort data=seds;
	by fipsst year;

*PUC data;

data puc (drop = fipsstn postalst);
	length FIPSST $2.;
	set exp.puc (keep = postalst year appoint elect); *Public Utility Commission Selection Data - see Table W1 in web appendix for details;
	rename appoint = puc_appoint;
	rename elect = puc_elect;

	fipsstn = stfips(postalst);
        FIPSST = translate(right(put(fipsstn, 2.)), '0', ' ');

	*Make correction to TN (see PUC Elect Variable Compare.doc);
	if postalst = 'TN' and year ge 1996 then do;
		appoint = 1;
		elect = 0;
	end;


proc sort data=puc;
	by fipsst year;

*NARUC data;

data naruc (drop = fips censt division region postalst state);
	set elec.naruc; *NARUC data - see pp. 5-6 of web appendix for details;

proc sort data=naruc;
	by fipsst year;

*Wide NARUC variables;

%w(1977)
%w(1982)
%w(1987)
%w(1992)
%w(1996)

data wide;
	merge w1977 w1982 w1987 w1992 w1996;
	by fipsst;

*Merge;

data pqem;
	merge pqem wide;
	by fipsst;

data styr_elas;
	merge pqem (in=a) seds (in=b) puc (in=c) naruc (in=d);
	by fipsst year;
	if a; *Keep only state-years in PQEM;

	label mean_elas_nf_pe = "Mean Elasticity of Electricity Price w.r.t. Quantity (Purchase-Weighted, No Utility FE)";
	label mean_elas_pe = "Mean Elasticity of Electricity Price w.r.t. Quantity (Purchase-Weighted, Utility FE)";
	label mean_elas_nf_tvs = "Mean Elasticity of Electricity Price w.r.t. Quantity (Shipments-Weighted, No Utility FE)";
	label mean_elas_tvs = "Mean Elasticity of Electricity Price w.r.t. Quantity (Shipments-Weighted, Utility FE)";

proc freq data=styr_elas;
	tables fipsst year;

proc print data=styr_elas (obs=5);
	title1 "Print 5 obs styr_elas";

*Fill in fac_i_adj and rank_s_adj to have values in every year;

data styr_elas;
	set styr_elas;
	*Set early years to zero;
	if year = 1963 or year = 1967 or year = 1972 then do;
		fac_i_adj = 0;
		rank_s_adj = 0;
	end;
	*Interpolate middle years;
	else if year ge 1973 and year le 1977 then do;
		yrd1 = year - 1972;
		facd1 = (fac_i_adj_1977 - 0) / 5;
		fac_i_adj = 0 + (facd1 * yrd1);
		rd1 = (rank_s_adj_1977 - 0) / 5;
		rank_s_adj = 0 + (rd1 * yrd1);
	end;
	else if year ge 1977 and year le 1982 then do;
		yrd2 = year - 1977;
		facd2 = (fac_i_adj_1982 - fac_i_adj_1977) / 5;
		fac_i_adj = fac_i_adj_1977 + (facd2 * yrd2);
		rd2 = (rank_s_adj_1982 - rank_s_adj_1977) / 5;
		rank_s_adj = rank_s_adj_1977 + (rd2 * yrd2);
	end;
	else if year ge 1982 and year le 1987 then do;
		yrd3 = year - 1982;
		facd3 = (fac_i_adj_1987 - fac_i_adj_1982) / 5;
		fac_i_adj = fac_i_adj_1982 + (facd3 * yrd3);
		rd3 = (rank_s_adj_1987 - rank_s_adj_1982) / 5;
		rank_s_adj = rank_s_adj_1982 + (rd3 * yrd3);
	end;
	else if year ge 1987 and year le 1992 then do;
		yrd4 = year - 1987;
		facd4 = (fac_i_adj_1992 - fac_i_adj_1987) / 5;
		fac_i_adj = fac_i_adj_1987 + (facd4 * yrd4);
		rd4 = (rank_s_adj_1992 - rank_s_adj_1987) / 5;
		rank_s_adj = rank_s_adj_1987 + (rd4 * yrd4);
	end;
	else if year ge 1992 and year le 1996 then do;
		yrd5 = year - 1992;
		facd5 = (fac_i_adj_1996 - fac_i_adj_1992) / 4;
		fac_i_adj = fac_i_adj_1992 + (facd5 * yrd5);
		rd5 = (rank_s_adj_1996 - rank_s_adj_1992) / 5;
		rank_s_adj = rank_s_adj_1992 + (rd5 * yrd5);
	end;
	*Set later years to their 1996 value;
	else if year > 1996 then do;
		fac_i_adj = fac_i_adj_1996;
		rank_s_adj = rank_s_adj_1996;
	end;

*Check;

proc print data=styr_elas (obs=100);
	var fipsst year fac_i_adj rank_s_adj dpf1;
	title1 "Check interpolation fac_i_adj and rank_s_adj";

**Electricity Intensity Measure;

*Get PQEM plant-level data;

data pqem2;
        set elec.gs6300 (keep = flag5_pw ppn year pe ee cp cr cw cf sw newind fipsst state wt tvs);

        rename state = postalst;

        if flag5_pw ne 1;

        *Create variable cost;
        vc = sum(cp, cr, cw, ee, cf, sw);

*Fipsst;

data state2;
        set pqem2 (keep = postalst fipsst);

proc sort data=state2 nodupkey;
        by postalst;

*Calculate the electricity cost share by 4-digit industry;

proc sort data=pqem2;
        by year newind;

proc summary data=pqem2;
        by year newind;
        output out=indyr
        sum(ee) = sum_ee
        sum(vc) = sum_vc
        ;

data indyr;
        set indyr;
        elec_share = sum_ee / sum_vc;

*Merge back to plant-level data;

data pqem2;
        merge pqem2 (in=a) indyr (in=b);
        by year newind;
        wttvs = wt * tvs;

*Create Index 1: Shipments-weighted mean by state of industry cost shares;

proc sort data=pqem2;
        by year postalst;

proc summary data=pqem2 vardef=wdf;
        weight wttvs;
        by year postalst;
        output out=ind1
        mean(elec_share) = index1
        ;

data ind1;
        set ind1 (keep = year postalst index1);

*Put sum_pe and sum_tvs on ind1;

proc sort data=ind1;
	by year postalst;

data temp_ind1;
	set styr_elas (keep = year postalst sum_pe sum_tvs);

proc sort data=temp_ind1;
	by year postalst;

data ind1;
	merge ind1 temp_ind1;
	by year postalst;

**Create demeaned version of index 1;

*Demeaned by year;

*Create annual means;

*Unweighted;

proc summary data=ind1;
	by year;
	output out=ind1yrmn
	mean(index1) = mean_yr_index1
	;

data ind1yrmn;
	set ind1yrmn (keep = year mean_yr_index1);

*Purchase-weighted;

proc summary data=ind1;
	by year;
	weight sum_pe;
	output out=ind1yrmn_pe
	mean(index1) = mean_yr_index1_pe
	;

data ind1yrmn_pe;
	set ind1yrmn_pe (keep = year mean_yr_index1_pe);

*Shipments-weighted;

proc summary data=ind1;
	by year;
	weight sum_tvs;
	output out=ind1yrmn_tvs
	mean(index1) = mean_yr_index1_tvs
	;

data ind1yrmn_tvs;
	set ind1yrmn_tvs (keep = year mean_yr_index1_tvs);

*Merge to ind1;

data ind1;
	merge ind1 ind1yrmn ind1yrmn_pe ind1yrmn_tvs;
	by year;

*Demeaned by grand mean;

*Create grand mean;

*Unweighted;

proc summary data=ind1;
	output out=ind1mn
	mean(index1) = mean_grand_index1
	;

*Merge to ind1;

data ind1 (keep = year postalst index1 index1_g index1_yr index1_yr_pe index1_yr_tvs);
	if _N_ = 1 then set ind1mn;
	set ind1;

	*Demean by grand mean;
	index1_g = index1 - mean_grand_index1;
	*Demean by annual mean;
	index1_yr = (index1 - mean_yr_index1) * 100;
	index1_yr_pe = (index1 - mean_yr_index1_pe) * 100;
	index1_yr_tvs = (index1 - mean_yr_index1_tvs) * 100;

*Merge to styr_elas;

proc sort data=styr_elas;
	by year postalst;

data styr_elas;
	merge styr_elas ind1;
	by year postalst;

**Fix caps for later ordering in csv;

data styr_elas;
	set styr_elas;
	rename PF1 = pf1t;
	rename PF2 = pf2t;
	rename PF3 = pf3t;
	rename SHARE_GEN_HY = share_gen_hyt;
	rename SHARE_GEN_NU = share_gen_nut;
	rename SHARE_GEN_OT = share_gen_ott;
	rename SHARE_GEN_CL = share_gen_clt;
	rename SHARE_GEN_PN = share_gen_pnt;
	rename SHARE_GEN_GE = share_gen_get;
	rename SHARE_GEN_NG = share_gen_ngt;
	rename SHARE_GEN_PA = share_gen_pat;
	rename SHARE_GEN_WW = share_gen_wwt;
	rename SHARE_GEN_WY = share_gen_wyt;
	

*Final styr_elas dataset;

data styr_elas;
	set styr_elas (drop = rd: yrd: facd: ddb_adj_1992 ddb_adj_1996 efr_adj_1992 efr_adj_1996 rank_s_short_adj_1992 
		rank_s_short_adj_1996);
	rename PF1t = pf1;
	rename PF2t = pf2;
	rename PF3t = pf3;
	rename SHARE_GEN_HYt = share_gen_hy;
	rename SHARE_GEN_NUt = share_gen_nu;
	rename SHARE_GEN_OTt = share_gen_ot;
	rename SHARE_GEN_CLt = share_gen_cl;
	rename SHARE_GEN_PNt = share_gen_pn;
	rename SHARE_GEN_GEt = share_gen_ge;
	rename SHARE_GEN_NGt = share_gen_ng;
	rename SHARE_GEN_PAt = share_gen_pa;
	rename SHARE_GEN_WWt = share_gen_ww;
	rename SHARE_GEN_WYt = share_gen_wy;

	label dpf1 = "One Year Log Change in Fossil Fuel Price Index pf1";
	label dpf2 = "One Year Log Change in Fossil Fuel Price Index pf2";
	label dpf3 = "One Year Log Change in Fossil Fuel Price Index pf3";

	label ddb_adj_1977 = "Adjusted: 1977 Discouraging Declining Block Rates (Yes = 1, No = 0)";
	label ddb_adj_1982 = "Adjusted: 1982 Discouraging Declining Block Rates (Yes = 1, No = 0)";
	label ddb_adj_1987 = "Adjusted: 1987 Discouraging Declining Block Rates (Yes = 1, No = 0)";

	label db_adj_1977 = "Adjusted: 1977 Declining Block Rates Approved (Yes = 0, No = 1)";
	label db_adj_1982 = "Adjusted: 1982 Declining Block Rates Approved (Yes = 0, No = 1)";
	label db_adj_1987 = "Adjusted: 1987 Declining Block Rates Approved (Yes = 0, No = 1)";
	label db_adj_1992 = "Adjusted: 1992 Declining Block Rates Approved (Yes = 0, No = 1)";
	label db_adj_1996 = "Adjusted: 1996 Declining Block Rates Approved (Yes = 0, No = 1)";

	label efr_adj_1977 = "Adjusted: 1977 Encouraging Flat Rates (Yes = 1, No = 0)";
	label efr_adj_1982 = "Adjusted: 1982 Encouraging Flat Rates (Yes = 1, No = 0)";
	label efr_adj_1987 = "Adjusted: 1987 Encouraging Flat Rates (Yes = 1, No = 0)";

	label fac_i_adj_1977 = "Adjusted: 1977 Fuel Adjsutment Clause - Interpreted (1 = Totally automatic, 0.75 = Partially 
		automatic, 0.5 = Allowed - not automatic, 0 = Not allowed)";
	label fac_i_adj_1982 = "Adjusted: 1982 Fuel Adjsutment Clause - Interpreted (1 = Totally automatic, 0.75 = Partially 
		automatic, 0.5 = Allowed - not automatic, 0 = Not allowed)";
	label fac_i_adj_1987 = "Adjusted: 1987 Fuel Adjsutment Clause - Interpreted (1 = Totally automatic, 0.75 = Partially 
		automatic, 0.5 = Allowed - not automatic, 0 = Not allowed)";
	label fac_i_adj_1992 = "Adjusted: 1992 Fuel Adjsutment Clause - Interpreted (1 = Totally automatic, 0.75 = Partially 
		automatic, 0.5 = Allowed - not automatic, 0 = Not allowed)";
	label fac_i_adj_1996 = "Adjusted: 1996 Fuel Adjsutment Clause - Interpreted (1 = Totally automatic, 0.75 = Partially 
		automatic, 0.5 = Allowed - not automatic, 0 = Not allowed)";

	label fr_adj_1977 = "Adjusted: 1977 Flat Rates Approved (Yes = 1, No = 0)";
	label fr_adj_1982 = "Adjusted: 1982 Flat Rates Approved (Yes = 1, No = 0)";
	label fr_adj_1987 = "Adjusted: 1987 Flat Rates Approved (Yes = 1, No = 0)";
	label fr_adj_1992 = "Adjusted: 1992 Flat Rates Approved (Yes = 1, No = 0)";
	label fr_adj_1996 = "Adjusted: 1996 Flat Rates Approved (Yes = 1, No = 0)";

	label inter_adj_1977 = "Adjusted: 1977 Interuptible Rates Approved (Yes = 1, No = 0)";
	label inter_adj_1982 = "Adjusted: 1982 Interuptible Rates Approved (Yes = 1, No = 0)";
	label inter_adj_1987 = "Adjusted: 1987 Interuptible Rates Approved (Yes = 1, No = 0)";
	label inter_adj_1992 = "Adjusted: 1992 Interuptible Rates Approved (Yes = 1, No = 0)";
	label inter_adj_1996 = "Adjusted: 1996 Interuptible Rates Approved (Yes = 1, No = 0)";

	label rank_s_adj_1977 = "Adjusted: 1977 Raw sum of db, inter, seas, fr";
	label rank_s_adj_1982 = "Adjusted: 1982 Raw sum of db, inter, seas, fr";
	label rank_s_adj_1987 = "Adjusted: 1987 Raw sum of db, inter, seas, fr";
	label rank_s_adj_1992 = "Adjusted: 1992 Raw sum of db, inter, seas, fr";
	label rank_s_adj_1996 = "Adjusted: 1996 Raw sum of db, inter, seas, fr";

	label rank_s_short_adj_1977 = "Adjusted: 1977 Raw sum of ddb, efr";
	label rank_s_short_adj_1982 = "Adjusted: 1982 Raw sum of ddb, efr";
	label rank_s_short_adj_1987 = "Adjusted: 1987 Raw sum of ddb, efr";

	label seas_adj_1977 = "Adjusted: 1977 Seasonal rates approved (Yes = 1, No = 0)";
	label seas_adj_1982 = "Adjusted: 1982 Seasonal rates approved (Yes = 1, No = 0)";
	label seas_adj_1987 = "Adjusted: 1987 Seasonal rates approved (Yes = 1, No = 0)";
	label seas_adj_1992 = "Adjusted: 1992 Seasonal rates approved (Yes = 1, No = 0)";
	label seas_adj_1996 = "Adjusted: 1996 Seasonal rates approved (Yes = 1, No = 0)";

	label sum_pe = "State-Year Sample-weighted Sum Quantity of Purchased Electricity";
	label sum_tvs = "State-Year Sample-weighted Sum Total Value of Shipments";

	label index1 = "State-Year Shipments-weighted Mean of National Industry Electricity Cost Shares";
	label index1_g = "State-Year Shipments-weighted Mean of National Industry Electricity Cost Shares 
(Demeaned by grand mean)";
	label index1_yr_pe = "State-Year Shipments-weighted Mean of National Industry Electricity Cost Shares 
(Demeaned by purchase-weighted annual mean) * 100";
	label index1_yr_tvs = "State-Year Shipments-weighted Mean of National Industry Electricity Cost Shares 
(Demeaned by shipments-weighted annual mean) * 100";
	label index1_yr = "State-Year Shipments-weighted Mean of National Industry Electricity Cost Shares 
(Demeaned by annual mean) * 100";


data e;
	set styr_elas;
	wtpe100 = sum_pe * 100;
	wttvs100 = sum_tvs * 100;
	
%d(1967)
%d(1982)

*Sum up to US total PE by year;

proc sort data=e;
	by year;

proc summary data=e;	
	by year;
	output out=ustot
	sum(sum_pe) = us_pe;

*Grab total for a given year;

%g(1967)

data eboth;
	merge e1967 e1982;
	by fipsst;
	delas_nf_pe = elas_nf_pe_1982 - elas_nf_pe_1967;
	delas_pe = elas_pe_1982 - elas_pe_1967;
	delas_nf_tvs = elas_nf_tvs_1982 - elas_nf_tvs_1967;
	delas_tvs = elas_tvs_1982 - elas_tvs_1967;

data eboth;
	if _N_ = 1 then set us1967;
	set eboth;

%r(elas_nf_pe)
%r(elas_pe)
%r(elas_nf_tvs)
%r(elas_tvs)
	
data e;
	merge eboth_elas_nf_pe eboth_elas_pe eboth_elas_nf_tvs eboth_elas_tvs;
	by postalst;

data rp001;
	set e;

*Get plant-level elasticities;

data outurall;
	set elec.outurall;
	pe = exp(logpe);

proc sort data=outurall;
	by fipsst;

data outurall_nf;
	set elec.outurall_nf;

proc sort data=outurall_nf;
	by fipsst;

*Get state ranks for top and bottom;

data rank;
	set rp001;

proc sort data=rank;
	by fipsst;

*Merge;

data outurall;
	merge outurall (in=a) rank (in=b);
	by fipsst;
	if a;

data outurall_nf;
	merge outurall_nf (in=a) rank (in=b);
	by fipsst;
	if a;

*Calculate weighted elasticities;

%rd1(outurall, tvs, elas)
%rd1(outurall, pe, elas)
%rd1(outurall_nf, tvs, elas_nf)
%rd1(outurall_nf, pe, elas_nf)

data rp002;
	merge felas_tvs felas_pe felas_nf_tvs felas_nf_pe;
	by year;

*Data for 2 top panels of Figure 7;

proc export data=rp002 outfile="csv/figure7_1a.csv"
	dbms=csv replace;

run;

*Data for use in creating Figure 8;

data elec.styr_elas;
	set styr_elas;

proc datasets lib=work kill;
run;
quit;

**Bottom panels of Figure 7;

*Get PQEM plant-level data;

data pqem;
	set elec.gs6300 (keep = flag5_pw ppn year pe ee cp cr cw cf sw newind fipsst state wt tvs);

	rename state = postalst;

	if flag5_pw ne 1;

	*Create variable cost;
	vc = sum(cp, cr, cw, ee, cf, sw);

*Fipsst;

data state;
	set pqem (keep = postalst fipsst);

proc sort data=state nodupkey;
	by postalst;
	

*Calculate the electricity cost share by 4-digit industry;

proc sort data=pqem;
	by year newind;

proc summary data=pqem;
	by year newind;
	output out=indyr
	sum(ee) = sum_ee
	sum(vc) = sum_vc
	;

data indyr;
	set indyr;
	elec_share = sum_ee / sum_vc;

*Merge back to plant-level data;

data pqem;
	merge pqem (in=a) indyr (in=b);
	by year newind;

*For this exercise we are interested in 1967 so keep only 1967;

data p1967;
	set pqem;
	if year = 1967;
	wttvs = wt * tvs;

*Create Index 1: Shipments-weighted mean by state of industry cost shares;

proc sort data=p1967;
	by postalst;

proc summary data=p1967 vardef=wdf;
	weight wttvs;
	by postalst;
	output out=ind1
	mean(elec_share) = index1
	;

data all;
	set ind1 (keep = postalst index1);

*Sort states into quartiles by index1;

%s(1)

*Create final dataset;

proc sort data=all;
	by postalst;

data all;
	merge all state;
	by postalst;

data rp008 (keep = fipsst postalst index: top: bottom: q1 fipsst);
	set all;

proc print data=rp008;
	title1 "Check final rp008";

*Get plant-level elasticities;

data outurall;
	set elec.outurall;
	pe = exp(logpe);

proc sort data=outurall;
	by fipsst;

data outurall_nf;
	set elec.outurall_nf;

proc sort data=outurall_nf;
	by fipsst;

*Get state ranks for top and bottom - ranks for each of three indices;

data rank;
	set rp008;

proc sort data=rank;
	by fipsst;

*Merge;

data outurall;
	merge outurall (in=a) rank (in=b);
	by fipsst;
	if a;

data outurall_nf;
	merge outurall_nf (in=a) rank (in=b);
	by fipsst;
	if a;

*Calculate weighted elasticities;

%d1_2(outurall, tvs, elas)
%d1_2(outurall, pe, elas)
%d1_2(outurall_nf, tvs, elas_nf)
%d1_2(outurall_nf, pe, elas_nf)

data rp009 (drop = _type_ _freq_);
	merge 
		felas_tvs1 felas_pe1 felas_nf_tvs1 felas_nf_pe1
	;
	by year;

proc export data=rp009 outfile="csv/figure7_2.csv"
	dbms=csv replace;

run;


